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Abstract. The theory of viscosity solutions has been effective for represent- 
ing and approximating weak solutions to fully nonlinear Partial Differential 
Equations (PDEs) such as the elliptic Monge-Ampere equation. The approx- 
imation theory of Barles-Souganidis [BS91] requires that numerical schemes 
be monotone (or elliptic in the sense of [Obe06]). But such schemes have 
limited accuracy. In this article, we establish a convergence result for filtered 
schemes, which are nearly monotone. This allows us to construct finite differ- 
ence discretizations of arbitrarily high-order. We demonstrate that the higher 
accuracy is achieved when solutions are sufficiently smooth. In addition, the 
filtered scheme provides a natural detection principle for singularities. We 
employ this framework to construct a formally second-order scheme for the 
Monge-Ampere equation and present computational results on smooth and 
singular solutions. 
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1. Introduction 

The numerical approximation of the Monge- Ampere equation is a problem of 
current interest, because of the many applications of the equation to various fields, 
and because the equation is the prototypical fully nonlinear elliptic Partial Differ- 
ential Equation (PDE). Thus building effective (convergent, fast, accurate) solvers 
to this equation demonstrates the possibility of effectively solving a wide class of 
fully nonlinear PDEs, which until recently were believed to be intractable. 

We consider the Monge- Ampere equation in a convex bounded subset £7 C M'' 

(MA) det{D^u{x)) = f{x), for x in fl, 

where det{D^u), is the determinant of the Hessian of the function u. We include 
Dirichlet boundary conditions on the boundary dfl, 

(D) u{x) = g{x), for x on dil. 

This equation is augmented by the convexity constraint 

(C) u is convex, 

which is necessary for the equation to be elliptic. 

This article builds on a series of papers which have developed solution methods 
for the Monge-Ampere equation. The foundation of schemes for a wide class of 
nonlinear parabolic and elliptic equations was developed in [Obc06]. The first 
convergent scheme for the Monge-Ampere equation was built in [ObeOSb]; this 
was restricted to two dimensions and to a slow iterative solver. Implicit solution 
methods were first developed in [BFOlO], where it was demonstrated that the use 
of non-monotone schemes led to slow solvers for singular solutions. In [FOlla] a 
higher dimensional monotone discretization was constructed, a fast Newton solver 
was also implemented. The convergent discretizations use a wide stencil scheme, 
which leads to a reduction on accuracy which reflects the directional resolution 
of the stencil. While this cannot be avoided on singular solutions, it is desirable 
to have a more accurate solver on (rare) smooth solutions. In [FOllb] a hybrid 
solver was built, which combined the advantages of accuracy in smooth regions, 
and robustness (convergence and stability) near singularities. However, this was 
accomplished at the expense of a convergence proof. While the Dirichlet problem is 
a natural starting point, for applications related to mapping problems or Optimal 
Transportation, other boundary conditions are used. These boundary conditions 
were implemented in [Frol2]. In a work in progress, the filtered scheme has also 
been apphed to the Optimal Transportation boundary conditions [BF012]. 
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1.1. Contribution of this work. In this article, we improve on previous results 
by building a convergent, higher order accurate scheme for the Monge- Ampere equa- 
tion, combining the advantages of the hybrid scheme in [FOllb] with the conver- 
gence proof in [FOlla]. Our result requires that we extend the convergence theory 
of Barles and Souganidis by considering the more general class of nearly monotone 
schemes. This extension is of independent interest, and applies to schemes for a 
wide class of elliptic equations. Once this theoretical result is established, it leads 
to a natural and simple method for constructing accurate discretizations of the 
Monge- Ampere equation, and indeed for the entire class of nonlinear elliptic PDEs, 
given the foundation of a monotone elliptic scheme. 

The combined schemes are called filtered finite difference approximations. We 
provide a proof that solutions of the filtered scheme exist and converge to the 
viscosity solution of the underlying PDE. The theory ensures, and computational 
results verify, that solutions of this scheme converge to the viscosity solution of the 
Monge- Ampere equation even in the singular setting. Using Newton's method, we 
construct a fast solver for the resulting nonlinear system of equations. 

The advantage of the more accurate scheme is obvious when solutions are smooth. 
In the case of a singular solution, the added accuracy of the filtered scheme is 
redundant, at least near singular parts of the solution. However, the since the low 
accuracy of the monotone scheme is not a problem on singular solutions, the filtered 
scheme allows for the use of a narrow stencil in general, while still achieving full 
accuracy on smooth solutions. 

1.2. The heuristics of a hybrid scheme. A natural way to build hybrid schemes 
is to make a convex combination of a stable, convergent scheme, Fm , and an accu- 
rate, less stable scheme. Fa- The weighting of the convex combination is determined 
by a function which measures the regularity of the solution so that 



where < < 1 is a continuous function which is 1 in a neighbourhood of a 
singularity and goes to elsewhere. In general Wg could be determined by checking 
the size of a derivative of u, for example ||£)^u|| > 1/h, where h is the grid spacing. 
We record the idea of a generic hybrid scheme with the following schematic 



So the function Wg depends on the solution, and this could potentially lead to 
instabilities unrelated to the stability of the underlying schemes Fm and Fa- 

In [FOllb], Fm corresponded to a wide stencil, monotone elliptic scheme, and Fa 
corresponded to the standard nine point finite difference scheme. We review those 
schemes below. In that article, we were able to use the regularity theory for the 
Dirichlet problem for (MA) to determine Ws a priori. However, in other problems 
of interest, for example with Optimal Transportation boundary conditions, Wg will 
depend on the solution. Since the resulting hybrid scheme was not monotone, we 
could not apply the Barles-Souganidis theorem and there was no other obvious way 
to prove convergence. However the method worked well in practice. 

To summarize: hybrid schemes are practical tools which blend the accuracy and 
stability of the underlying schemes. But they are defined in an ad hoc manner, and 
there is no guarantee of the observed stability and accuracy. 



Fh = WsFm + (1 - Ws)Fa, 



(1.1) 




Fm near sin gularity, e.g. \\D'^u\\ > 1/h, 
Fa elsewhere . 
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1.3. A motivating example for a filtered scheme. In this section, we present 
a filtered scheme on a simpler equation, to present the ideas which follow more 
clearly. 

The filtered scheme provided an intrinsic method for defining based on the size 
of the difference between the monotone operator and the accurate operator. For 
illustration purposes, consider the model equation 

F[u]{x) = - 1 

on the domain [—1,1] with boundary conditions u(— 1) = w(l) = 1. Then the 
viscosity solution is simply u{x) = \x\. Define the monotone upwind scheme 

which is first order accurate: F^{(j)) — F{(f)) = 0{h), for smooth 0. The scheme 
F^^ is consistent with the method of characteristics, and it satisfies a maximum 
principle, so solutions converge. The equation has the explicit form 

u{x) — max{u(a; + h), u{x — h)} — h. 
Next define the second order accurate, but unstable, centred difference scheme 

r.h \u{x + h) - u{x - h)\ 

- 2h ^' 

F^(0) - F{cj)) = 0{h^), for smooth </>. 

A natural definition of a singularity of the equation is when \ux:x\ is large. For 
the finite difference scheme, this can be interpreted as 

\u{x + h)- 2u{x) + u{x-h)\ 1 
J? - h' 

Taking the difference between the two schemes, we obtain a similar condition, up 
to the scaling in h, 

ijph p/i I 1 \u{x + h) - 2u{x) + u{x - h)\ 
\Fa - Fm\ - 2 ^ • 

So a singularity can be defined by the condition 

1 



(1.2) 



F^ - F^, 



> . 

- h 



The inequality above leads to a hard condition. Instead we replace it with a soft 
condition using a continuous filter function. 

Definition 1 (Filter function). We define a filter function to be a continuous, 
bounded function, 5, which is equal to the identity in a neighbourhood of the 
origin and zero outside. 

For example, 
(1.3) Six) = 

See Figure 1. 
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Figure 1 . Filter function 



Then the filtered scheme can be defined as follows 



Fs = Fm + hS 



( 



Fa - Fm 
h 



) 



By virtue of the condition (1.2), the filtered scheme will be consistent with the 
condition for a hybrid scheme (1.1). (In general, for (1.2) to hold, the scaling h 
may depend on the order of accuracy of the schemes Fa and Fm)- 
By the definition (1.3), 



Thus the difference between the filtered scheme, Fg, and the monotone scheme, 
Fm, goes to zero with h, uniformly over all functions. This property is what allows 
us to prove convergence, since the filtered scheme is nearly monotone. 

1.4. Introduction to numerical methods for degenerate elliptic PDEs. 

There are two major challenges in building numerical solvers for nonlinear and 
degenerate elliptic PDEs. The first challenge is to build convergent approxima- 
tions, often with finite difference schemes. The second challenge is to build efficient 
solvers. 

The approximation theory developed by Barles and Souganidis [BS91] provides 
criteria for the convergence of approximation schemes: schemes that are consistent, 
monotone, and stable converge to the unique viscosity solution of a degenerate 
elliptic equation. However, this work does not indicate how to build such schemes, 
or how to produce fast solvers for the schemes. It is not obvious how to ensure that 
schemes satisfy the required comparison principle. The class of schemes for which 
this property holds was identified in [Obe06], and were called degenerate elliptic, 
by analogy with the structure condition for the PDE. 

An important distinction for this class of equations is between first order (Hamilton- 
Jacobi) equations and the second order (nonlinear elliptic) case. The theory of vis- 
cosity solutions [CIL92] covers both cases, but the corresponding numerical methods 
are quite different. In the first order case, where the method of characteristics is 
available, there are formulas for exact solutions (e.g. Hopf-Lax) and there is a con- 
nection with one-dimensional conservation laws [Eva98] . The second order case has 
more in common with divergence-structure elliptic equations; however, for degen- 
erate of fully nonlinear equations, many of the tools from the divergence-structure 
case (e.g. finite elements, multi grid solvers) have not been successfully applied. 



\Fs-Fm\ < 2h. 
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Much of the progress on discretization techniques and fast solvers is hmited to 
the first order case. For Hamilton- Jacobi equations, which are first order nonlinear 
PDEs, monotonicity is necessary for convergence. Early numerical papers stud- 
ied explicit schemes for time-dependent equations on uniform grids [CL84, Sou85]. 
These schemes have been extended to higher accuracy schemes, which include sec- 
ond order convergent methods — the central schemes [LTOO] , as well as higher order 
interpolation methods — the ENO schemes [0S91]. Semi-Lagrangian schemes take 
advantage of the method of characteristics to prove convergence [FF02] . These have 
been extended to the case of differential games [BFS99]. Two classes of fast solvers 
have been developed: fast marching [Set99, Tsi95] and fast sweeping [TCOZ03]. 
The fast marching and fast sweeping methods give fast solution methods for first 
order equations by taking advantage of the method of characteristics, which is not 
available in the second order case. 

There is much less work available for second order degenerate elliptic equa- 
tions. For uniformly elliptic PDEs, monotone schemes are not always necessary 
for convergence (for example, most higher order finite element methods are not 
monotone). However, for fully nonlinear or degenerate elliptic equations, the only 
convergence proof currently available requires that schemes be monotone. One way 
to ensure monotonicity is to use wide stencil finite difference schemes; this has 
been done for the equation for motion by mean curvature [Obe04], for the Infinity 
Laplace equation [Obe05] , for functions of the eigenvalues [ObeOSb] , for Hamilton- 
Jacobi-Bellman equations [BZ03], and for the convex envelope [ObcOSa]. Even 
for linear elliptic equations, a wide stencil may be necessary to build a monotone 
scheme [MW53]. 

1.5. Higher order numerical methods. Monotone schemes have limited accu- 
racy. For first-order equations, a consistent, monotone scheme will be at most 
first-order accurate. For second-order equations, the accuracy is at most second- 
order [Obe06]. Moreover, as we have already noted, wide stencils are needed to 
build monotone schemes for some equations. In this case, the formal accuracy of 
the scheme depends not only on the spatial resolution of the scheme, but also on 
the angular resolution of the scheme. To make large computations practical, it is 
desirable to restrict schemes to a reasonably narrow stencil width. However, this 
can place severe restrictions on the accuracy that can be achieved. 

For first-order Hamilton- Jacobi equations, some progress has been made in the 
construction of convergent higher-order schemes. Lions and Souganidis [LS95] con- 
sidered using a general higher-order scheme, which must be "filtered" to ensure the 
preservation of fractional one-sided second derivative bounds. These bounds, which 
are used instead of a monotonicity condition, coincide with the condition needed 
to ensure uniqueness of almost everywhere solutions of Hamilton-Jacobi equations; 
this approach does not generalize to second-order equations. Using an approach 
more similar to the one presented in this article, Abgrall [Abg09] proposed a "blend- 
ing" of a monotone and a higher-order scheme. Provided a solution to the blended 
scheme exists, it converges to the viscosity solution of the steady Hamilton-Jacobi 
equation. 

In the case of second-order fully nonlinear elliptic equations, we are not aware 
of any convergent higher-order schemes. 
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1.6. Numerical methods for the Monge-Ampere equation. For singular so- 
lutions, Newton's method combined with a non-monotone discretization of Monge- 
Ampere can become unstable [FOllb], which necessitates the use of other solvers 
for the discrete system of nonlinear equations. The simplest solver is to simply 
iterate the discrete version of the parabolic equation 

ut = det{D^u) - f 

using forward Euler. This method is slow, since the convergence time depends on 
the nonlinear CFL condition, which become more stringent with increasing prob- 
lem size [BFOlO]. Semi-implicit solvers can improve the solution time for smooth 
solutions, but can slow down or break down on singular solutions [BFOlO]. Similar 
behavior was observed using other discretizations [DG06]. 

However, in addition to suffering from a severe time-step restriction, this ap- 
proach does not enforce the convexity constraint. In particular, this type of itera- 
tion is unstable for solutions that are non-strictly convex. In two dimensions, we 
constructed solution methods that selected the convex solution. While no conver- 
gence proof was available, the methods appeared to converge to the correct weak 
solution of the equation. However, convergence was very slow for non-smooth or 
non-strictly convex solutions. 

In order to ensure convergence, more sophisticated techniques are needed. An 
early work by Oliker and Prussner [OP88] described a method that converges to 
the Aleksandrov solution in two dimensions. This early method was used to solve 
a problem with only about a dozen grid points. Several other methods have been 
proposed in recent years [BGNSll, DG06, FN09, LR05]; these are similar to our 
standard finite difference methods in terms of the lack of convergence theory and 
behaviour on singular solutions. We draw particular attention to the vanishing 
moment method proposed by Feng and Neilan [FN09], which involves regularizing 
the equation by adding a small multiple of the bilaplacian. This approach will be 
discussed in more detail in §3.3. 

1.7. The finite difference discretization of Monge-Ampere. We describe the 
elliptic (monotone) representation of the Monge-Ampere operator used in [FOlla, 
Frol2]. 

We begin with some informal remarks to explain the representation of the op- 
erator. It bears repeating that the standard finite difference discretization is not 
elliptic or convex. 

One difficulty comes from the off diagonal terms u^y. So by using wide stencils 
and rotating the coordinate system we can hope find the coordinate system in 
which the Hessian matrix is diagonal. But finding this coordinate system must be 
done in a monotone manner. The way to do this is to use a version of Hadamard's 
inequality, which we interpret to give an expression for the determinant of a positive 
definite matrix as a minimization. 

Hadamard's inequality states that for a positive semidefinite d x d matrix M , 

d 

det(M) < Y[mu 

i=l 
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We can force equality by choosing a coordinate system where M is diagonal. We 
write this as 



det(Af) = min J^nf Mrii, for M positive definite 

^'^'^ i=l 

where O is the set of all orthogonal d x d matrices, and ni,i ~ 1, . . . d are the rows 
of N. 

Another difficulty comes from the fact that the operator is elliptic only when 
M is positive definite. So we would like a continuous extension of the operator to 
non-positive definite matrices. A first possibility is 



d 

det'(M) = min TT max {nfMni,0\ 



4 = 1 



which extends the formula. However, we would also like to apply Newton's method, 
which computes (and inverts) the gradient of the operator. But this operator has 
a zero gradient. So instead, we use the following extension, which gives a negative 
value with a non-zero gradient on (most) non-positive definite matrices. 



d 



det(M) = min < Y[ ™ax {njMm, 0} + ^ min {nf Mn,, 0} 



^4 = 1 4=1 



Neo 
where 

d~~t(M) Met(Af), M positive semi-definite 

I Ai + • • • + Afc, Ai < • • • < Afc < are the negative eigenvalues of M 

We now apply the representation to the matrix D^u{x). The term nj D'^urii has a 
natural interpretation as a directional second derivatives, in the direction n^. This 
allows us to write 

f d^u \ . f d'^u 



det(-D'^it) — min < max < -r^, > 4- min < 
' Neo\ Xdnf J 1 



2 ' 



t dn, 

where O is the set of all orthogonal d x d matrices, and ni,i = 1, . . . d are the rows 
of N. 

The next stage is to replace the derivatives with finite differences, to arrive at a 
semi-discrete representation. However, we have only a limited number of directions 
available on the grid. So we limit ourselves to considering a finite number of vectors 
rii that lie on the grid and have a fixed maximum length; this is the directional 
discretization, which gives us the angular resolution dO of our stencil; see Figure 2. 
In this figure, values on the boundary are used to maintain the directional resolution 
dO (at the expense of lower order accuracy in space because the distances from the 
reference point are not equal). Another option is to use narrower stencils as the 
boundary is approached, which leads to lower angular resolution, but better spatial 
resolution. We denote the resulting set of orthogonal vectors by G- 

Each of the directional derivatives in the Monge- Ampere operator is then dis- 
cretized using centered differences: 

Vi.^Ui = }^ {u{xi + vh) + u{xi - vh) - 2u{xi)) . 
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Figure 2. Neighboring grid points used for width one (green), 
two (yellow), and three (blue) stencils. The illustration shows the 
neighbors in the first quadrant. The modification near the bound- 
ary is illustrated in the second and third figures. 



Thus the discrete version of the Monge- Ampere operator is 

{d d 1 

TT maxlV^^ u„ 6} + y^ minlV^,^^ u„ 6} } 
J=i U J 

where (5 > is a small parameter used to bound the second directional derivatives 
away from zero. 

2. Viscosity Solutions of Degenerate Elliptic Equations 

2.f . Background from viscosity solutions. In this section we review the theory 
of viscosity solutions. See the references above and [CIL92]. 

Given an open domain fl C M.'^, a second order partial differential equation on 

is a function 

which we write as F{x,r,p, M). Given a function u G C^(ri) we let p — Vu and 
M — D^u denote the gradient and Hessian of u, respectively. We include Dirichlet 
boundary conditions into the operator F in order to pose equation (PDE) in the 
closed domain by defining 

F[u]{x) — u{x) — g{x), X e d^l. 

Then the function u E C^{il) is a solution of the PDE in il if 

(PDE) F[u]{x) ^ F{x,u{x),\7u{x),D^u{x)) ^0, x e il. 

However, most of the PDEs we consider fail to have classical solutions under general 
assumptions on the data. This motivates the definition of viscosity solutions. 

Definition 2. The equation (PDE) is degenerate elliptic if 

F{x,r,p,X) <F{x,s,p,Y) 

for all X e 0., r,s e R, p e E", X,Y e with X >Y and r < s. Here X >Y 
means that X — Y is a, positive definite matrix. 
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The definition is related to the comparison principle. 

Lemma 1. The function F{x,r,p, X) is degenerate elliptic if and only if whenever 
a; is a non-negative local maximum of m — u, for u,v G C^, F[u]{x) > 

Proof. If x is a local maximum, u > v, Dv = Du, and D^u < D^v, at x. Then 
F{x, u, Du, D^u) = F{x, u, Dv, D^u) > F{x, v, Dv, D^u) > F{x, v, Dv, D^v). Here 
we have used the definition degenerate elliptic. □ 

Elliptic equations need not have smooth solutions, which necessitates some no- 
tion of a weak solution. We are interested in the viscosity solution [CIL92]. Before 
we define this, we introduce the upper and lower semi-continuous envelopes of a 
function. 

Definition 3 (Upper and Lower Semi-Continuous Envelopes). The upper and lower 
semi-continuous envelopes of a function u{x) are defined, respectively, by 

u* (x) = limsupu(t/), 
M* (x) — lim inf u(y) . 

Definition 4 (Viscosity Solution). An upper (lower) semi-continuous function u is 
a viscosity sub(super)-solution of (PDE) if for every S C^(i7), whenever u — (p 
has a local maximum (minimum) at x e f2, then 

F^{x,u{x),\7<j){x),D^(j){x)) < 

{F*{x,u{x),V(j){x),D^(l3{x)) > 0). 
A function u is a viscosity solution if it is both a sub- and a super-solution. See 
Figure 3. 




Figure 3. Touching u{x) by a smooth test function 



Remark 1 . When checking the definition of a viscosity solution, we can limit our- 
selves to considering unique, strict, global maxima (minima) oi u — 4> with a value 
of zero at the extremum. See, for example, [Koi04, Prop. 2.2]. For an accessible 
introduction to viscosity solution in the first order case, see [Eva98, Chapter 10]. 

The equations we consider satisfy a comparison principle. If u e USC(n) is a 
sub-solution and v G LSC(r2) is a super-solution of (PDE), then u < v on Cl. The 
proof of this result is one of the main technical arguments in the viscosity solutions 
theory [CIL92]. For a general proof and also simplified proofs in special cases, we 
refer to [Koi04, Chapter 3]. 
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2.2. Approximation schemes. An approximation scheme is family of functions 

parameterized by e G M+, 

which we write as F'^{x, r, Given a function v £ L°°{Q), we write 

(2.1) F'[v]ix)^F'{x,v{x),v{-)) 

where e > is small. The function u"^ is a solution of the scheme F'^ if 
(PDE)<^ = 0, for all x in 

Remark 2. For finite difference schemes, which are defined on a grid, piecewise linear 
interpolation is monotone, so we can extend the function continuously onto the 
domain. We assume in the current setting that we are working with the extended 
function. 

The important properties of a scheme follow. 

Definition 5 (Consistent). The scheme (2.1) is consistent with the equation (PDE) 
if for any smooth function (j) and x G Cl, 

limsup F^(y, <j>iy) + ^, 0(-) + C) < F*{x, c^{x) ,V c^{x) , c^[x)) , 
liminf F'{y, c^{y) + 0(-) + > i^. <t>ix), V c^[x) , cj>(x)) . 

Accuracy is defined inside the domain. To be more precise, we could also define 
accuracy at the boundary. 

Definition 6 (Accurate). The scheme (2.1) is accurate to 0{t°') if for any smooth 
function (j) and x £0,, 

F^[4>\-F[4>\^0{e"). 

Definition 7 (Stable). The scheme (2.1) is stable if any solution of (2.1) is 
bounded independently of e. 

Definition 8 (Monotone). The scheme (2.1) is monotone if for every e > 0, a; € H, 
s £ M and bounded u, u, 

u>v=^ F'(a;,s,u(-)) < i^'(a;,s, «(■))• 
Definition 9 (Elliptic). The scheme (2.1) is elliptic if it can be written 

F^[v\ = F%x,v{x),v{x)-v{-)), 
where F"^ is nondecreasing in its second and third arguments, 

(2.2) s < t, u{-) < «(•) =^ F%x,s,u{-)) < F'{x,t,v{-)) 
Elliptic schemes are monotone, since 

u(-) > u(-) s — u{-) < s — v{-) 

=^ F'{x, s,s~ u{-)) < F^{x, s,s- v{-)). 

In addition, it is shown in [Obc06] that under mild technical conditions solutions 
of elliptic schemes exist, satisfy a comparison principle, and are stable. These 
technical conditions are easily satisfied by finite difference schemes. For example if 
we consider the Dirichlet problem, a small perturbation of the scheme (which can 
go to zero with e) can be introduced to ensure that these conditions are satisfied. 
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A more complete theory for well-posedness of elliptic schemes can be found in 
[M0S12], which does not require perturbations. 

Definition 10. The scheme is a perturbation if there is a nonnegative modulus 
function m : M+ ^ M+ with 

lim m(e) — 

such that 

sup sup |F^[u](a;)| < m{e). 

Definition 11 (Nearly Monotone). The scheme F'^ is nearly monotone if it can be 
written as 

(2.3) F^[u]^FlAu]+F^^[u] 
where F^j is monotone and F^ is a perturbation. 

2.3. Convergence proof. Theorem 1, generalizes the corresponding result for 
monotone schemes in [BS91]. It uses Lemma 2 below, which is a standard re- 
sult from viscosity solutions theory. We include the proof since a reference to the 
result is not easily found. 

Theorem 1 (Convergence of Approximation Schemes). Let u be the unique solution 
of the degenerate elliptic (PDE). For each e > 0, let m*^ be a stable solution of the 
nearly elliptic approximation scheme (PDE)*^. Then 

u'^ — >■ u, locally uniformly, as e — >■ 

Proof of Theorem 1. Define 

u{x) = limsup u'(y) G USC{VL), 

u{x) = liminf u"(?/) G LSC{n). 

e— fO,y— 

These functions are bounded by the stability of solutions (Definition 7) . We record 
the fact that 

(2.4) u<u. 

If we know that u is a subsolution and u is a supersolution, then we could apply 
the comparison principle for (PDE) to u and u to conclude that 

u < u. 

Together these inequalities imply that u ^ u. The local uniform convergence follows 
from the definitions. 

It remains to show that m is a subsolution and u is a supersolution. 

We proceed to show that m is a sub-solution by testing the condition in defini- 
tion 4. Given a smooth test function let xq be a strict global maximum oi u — (j) 
with 4>{xq) — u{xq). 

By Lemma 2 (below), we can find sequences with 

< Vn 

where y„ is a global maximizer of u*^" — (j). 
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Define 

(2.5) ^« = ■«'"(yn) - 0(2/n)- 

Then 

£,n u{xo) - 0(a;o) = 
and u'^"{x) — 4>{x) < u'^"(?/„) — 4>{y.n) — £,n for any a; G In particular 

(2.6) ^,^.(.)_0(.)<^„. 

As a consequence of Definitions 8 and 11 we see tfiat 

(2.7) u{-)<vi-) F^x,s,u{-)) > F'{x,s,v{-)) -2m{e), 

for every e > 0, x G Cl, s E M.. 

Using the definitions above and (2.7), we find that 

= ](?/„) = i^'"(y„,u'" (?/„), u'"(-)), since u''" is a solution, 

= F^"(2/„,0(2;„)+en,0(-) + (^^'"(-)-'^(-))) by (2.5), 
> F'" (2/„, 0(?;„) + 0(-) + 60 - 2m{e) by (2.6) applied to (2.7). 

Next, we compute 

> liminf {^^^"(y„, ,/)(y„) + Cn, 0(-) + 6.) - 2m(e)} 
> hminf i^^(y,0(2;) + ^,0(.)+^) 

— F^{xct, (j){xo), \i'(f>{xn), D'^(f>{xn)) by definition 5 

= F^{xo,u{xo),V(f)ixQ),D'^(f){xo)), 

which shows that u is a subsolution. 

By a similar argument, we can show that u is a super-solution. □ 

Lemma 2 (Stability of Maxima). Suppose the family u'^ is bounded uniformly in e. 
Define 

u{x) = limsup u'(y) £ USC{n), 

Given a smooth function 0, let xq be a strict global maximum of u — t/i. Then there 
exist sequences: 

'e„ ^ 
2/n -J' xo 

^W'^iVn) u{xo) 
where y„ is a global maximizer of u*^" — 0. 

Proof of Lemma 2. From the definition of the limit superior, we can find sequences 

e„ -> 0, Zn ^ Xo 

such that 

Now we define j/„ €E to be maximizers of u'^" (x) — (j){x). 
We have 

u^^ivn) ~ (f>iyn) > u^'^izn) - 0(z„) '«(a;o) - 4'{xo) = 0. 
Also, for any 6 > and large enough n, 

(yn) - (t^iVn) < uiVn) " (/-(yn) + S < u{xo) ~ (/>(a;o) + (5 = (5. 
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Thus we have 

- HVu) 0. 

Now suppose we do not have j/„ — ^ a;o. Then (possibly through a subsequence) 
there is an i? > such that 

Also, since the maximum is strict, global, and unique, there is a > such that 

u{y) - Hy) <~K<o 

whenever |y — xo| > R. 

Thus for any 6 > and large enough n, 

{Vn) - (piVn) < u{yn) - 0(y„) + S < - K + S ^ - K < , 

which contradicts the fact that u'^"{yn) — <(>{yn) 0. We conclude that 

Finally, it is clear that 

|w'"(2/n) - u{xq)\ = |u'"(yn) - <P{xq)\ 

< l"'"(yn) - Hyn)\ + \Hyn) - 0(a;o)| 

^ 0. 

Therefore, 

u'^'^iyn) ~> u{xo)- □ 

2.4. Solutions of nearly monotone finite difference methods. Theorem 1 
assumed existence and stability of solutions to nearly monotone schemes. We show 
next that this follows from well-posedness of the underlying monotone schemes. 

For simplicity, we work with grid functions, i.e. we assume that is a finite 
set, for example a finite difference grid, which is identified with M^. For general 
approximation schemes, some other form of compactness can be used to achieve 
the same results. 

Lemma 3 (Existence and Stability of nearly monotone schemes). Suppose that 
solutions exist and are stable for the inhomogeneous problem for monotone scheme 

If is a continuous perturbation, then stable solutions exist for the nearly mono- 
tone scheme (2.3). 

Proof. Existence Fix e > 0. Write u — S{g) for the solution operator of the scheme 

By assumption, F^ is uniformly bounded on its domain, 

Fn{u) < R, for aU u. 
Since S is continuous, we have 

l|5(F^)ll < R2- 
In particular, for the ball Br^ C M^, 

S{F^^{Bn,)) c Bn,. 
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Applying Brouwer's fixed point theorem, we concluded that there is a fixed point 

S{Ff^(u*)) — u* , which means 

Thus we have identified a solution of the nearly monotone scheme. 

Stability We wish to show that any solution m*^ of the nearly monotone scheme 
can be bounded uniformly as e — 0. Then u'^ is a solution of 

regarding v as the unknown. Since is bounded, continuity of the solution 
operator for Flj gives a uniform bound for v. 

□ 



3. Applications 

3.1. Filtered schemes. We will now construct a more accurate nearly monotone 
scheme, starting from the monotone scheme F^j and the consistent (and more 
accurate) scheme F\. 

This is accomplished by setting 



(3.1) F^M=e"^ 



in (2. 3), where S{x) is a filter function, as defined in Definition 1. Since = 
e", (2.3) is a nearly monotone scheme. 

Lemma 4 (Filtered scheme is accurate). Suppose that the formal discretization 
errors of the schemes Fa, Fm are 0{e^*) and 0{e^'-') respectively. Choose the 
parameter a so that I3a> f^M > oi > Q. If is smooth, then — F%[(l)]. 

Proof. If (j) is smooth, then 

Fim-Flim ^ g(6^-) + 0(e^-O ^ Of,^--) < Oil) 

Using (3.1), we find that = F%[(j)\ - [</>], which means F''[(j)\ = F%[(j)\. □ 

Remark 3 (Singularity detection) . The filtered scheme provides an intrinsic singu- 
larity detection mechanism, which is adapted to the discretization (and the equa- 
tion) itself. This means fewer points are considered singular, compared to simply 
using a regularity condition. The selection principle is illustrated in Figure 6. 

3.2. Perturbations. Given the degenerate elliptic (PDE), and a continuous oper- 
ator G[u\, we can consider the nearly monotone approximation 

F'[u] = F[u]+eCS{G[u]/C), 

for some large constant C. Theorem 1 gives a convergence proof. We consider the 
special case of a solution u of F[v\, where G[u] is bounded by C. Then for e small 
enough we recover the unfiltered perturbation. 

For numerical methods, if we know a priori that the solution satisfies the bounds 
above, we can use an accurate method. However, for weak solutions of the PDE, 
this is not the case. 
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3.3. Bilaplacian regular izat ion. The Vanishing Moment Method proposed by 
Feng and Neilan [FN09] for solving the second-order eUiptic equation (PDE) where 
we specify the boundary conditions 

u = g, on 9ri 

by perturbing with the bilaplacian 

(3.2) F[u]+eA'^u^O. 

This also introduces additional boundary conditions, for example 

Au = 0, on dn. 

Here we provide a proof of convergence of the filtered PDE in the case of smooth 
solutions. 

Consider instead the filtered PDE 

(3.3) F'[u] = F[u]+eCS( — 
along with filtered boundary conditions 

(3.4) eS (^^^ = 0, on dQ, 

for a large constant C. The scheme we defined is nearly monotone. If | A^it| , \ Au\ < 
C then we recover (3.2). This follows since if the unperturbed equation is smooth, 
adding eA'^u only improves regularity in the interior. 

Theorem 1 then provides a convergence proof for the method in the case of 
smooth solutions. 

4. Computational Results 

In this section we implement the filtered scheme (3.1) for the Monge- Ampere 
equation, using a monotone scheme and standard finite differences. 

4.1. Numerical implementation. In this section, we implement a convergent 
filtered finite difference scheme for the Monge-Ampere equation. The purpose of 
this section is to demonstrate that the filtered scheme: 

(1) Allows for higher accuracy than the monotone scheme. 

(2) Can be solved efficiently using Newton's method. 

In practice, we do not expect the results to be much better than for the hybrid 
scheme. The improvement is that the hybrid scheme was overly conservative, so 
the filtered scheme reduces to the accurate scheme. But the real advantage of this 
method is that it can be used in situations (such as the Optimal Transportation 
problem) where the hybrid scheme may not converge. 

We now solve the Monge-Ampere equation using the filtered scheme 



e{h,d0) 



F^^'^^M = F^f[u] + e{h, d0)S 




Here Fm is the monotone scheme described in subsection 1.7 and Fa is the for- 
mally second order accurate scheme obtained using a standard centred-diS'erence 
discretisation (see [FOllb] for details). In addition, we have regularized the max 
operator to allow the gradient to be computed, as in [FOllb]. In the computations 
below, we used e{h, d0) ^Vh + d0/lO. 
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We solve the resulting system of nonlinear equations using Newton's method: 
u"+i = It" - {VFhK])-^FhK] 
where the Jacobian will be given by 

VFh[u] = (1 - S'[u]) VFm[u] + S'[u]VFa[u]. 
The derivative of the filter (1.3) is given by 

fl \x\<l 
S'{x) = 1 1 1< |x| < 2 
[o |a;|>2. 

However, allowing this derivative to take on negative values can lead to poorly 
conditioned or ill-posed linear systems. Instead, we approximate the Jacobian by 

VFh[u] = (1 - S"[m]) VFm[m] +max{S"[M],0}VFA[ii]. 

To properly assess the speed and accuracy of our filtered method, we also solve 
the Monge- Ampere equation using 

(1) The monotone scheme described in subsection 1.7. 

(2) The a priori hybrid scheme presented in [FOllb] (which has no convergence 
proof) . 

(3) The formally second-order standard scheme solved using one of the methods 
of [BFOlO]. 

4.2. Computational examples. We use our method to compute the following 
four representative examples described in [FOlla]. 

Throughout these definitions, we use x = (x, y) to denote a general point in 
and xq = (0.5, 0.5) for the center of the domain. 

Given a solution u(x, y), we obtain the mapping Vu(x, y) from the computational 
square to the image of the square. We visualize the optimal mappingFigure 4 by 
plotting the images of the constant x and constant y lines from the square on the 
image set. 

The first example, which is smooth and radial, is given by 

(4.1) «(x) = exp (^^^^^^ , /(x) = (l + ||x-xo|p)exp(||x-xo|p). 
The second example, which is C^, is given by 

(4.2) ^(x) = ^((||x-xo||-0.2)+)^ /(x)=(^l- 

The third example is smooth in the interior of the domain, but has an unbounded 
gradient near the boundary point (1, 1). The solution is given by 

(4.3) 7.(x) = -v/2-||x||2, /(x) = 2(2-||x||2)-^ 
The final example is the cone, which is only Lipschitz continuous. 

(4.4) u(x) ||x-xo||, / = /i = 7r(5x„. 

In fact, this solution is not even a viscosity solution; it must be understood using 
the more general notion of an Aleksandrov solution. In order to approximate the 
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solution on a grid with spatial resolution h using viscosity solutions, we approximate 
the measure /i by its average over the ball of radius h/2, which gives 



4//i2 for ||x-xo|| < /i/2, 
otherwise. 



The solutions and the corresponding mappings are plotted in Figure 4. 




(B) 

Figure 4. Top row: surface plot of the solutions, u{x, y). Bottow 
row: Mappings Wu{x, y) visualized by plotting the image of the 
constant x and y line segment in source domain. The solutions 
correspond to (from left to right) solution, solution, solution 
with blow-up, and cone solution (no mapping in this case, since it 
is singular). 



All computations are performed on the domain [0, 1] x [0, 1], which is discretized 
on an X grid. We choose e{h, d9) = y/h + d9 /lO. Computations are done using 
the 9, 17, and 33 point stencils. 

4.3. Discussion of computational results. Accuracy We begin by comparing 
the maximum error in the filtered scheme to the error in the monotone scheme. The 
accuracy is illustrated in Figure 5 as well as in Table 1. As we expect, the filtered 
scheme results in improved accuracy. The number of grid points was chosen so that 
in many examples the directional resolution error dominates the spatial resolution 
error. This is the reason that the accuracy of the monotone scheme tapers off. On 
the other hand, the hybrid and filtered schemes achieve better accuracy. For the 
cone example, the theory does not ensure convergence, and in this case, the standard 
scheme is more accurate, but the computation time for the standard scheme is very 
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long, and there are no guarantees for this scheme: the accuracy was much worse 
on other examples. 

Not surprisingly, the gains in accuracy are greatest on the smoothest solutions, 
when we can reasonably expect the more accurate scheme to be valid. On more 
singular examples, the standard scheme may not be valid, and the filtered scheme 
can be forced to choose the less accurate monotone scheme at many points in the 
domain. This is not a limitation of the scheme, though, since choosing the formally 
more accurate scheme everywhere can lead to instabilities or convergence to the 
wrong solution. 

We also compare our results to the a priori hybrid scheme, which, despite lacking 
a proof of convergence, was experimentally found to be efficient and more accurate 
than the monotone scheme [FOllb]. 
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Figure 5. Accuracy of the different schemes. Solutions are or- 
dered as in Figure 4. 



Figure 6 shows the weights of the monotone and accurate schemes determined 
by the filtered schemes. For the smooth solution, the accurate method is always 
used. Near singularities, the monotone method is selected. For comparison, the 
corresponding weights for the hybrid scheme are also shown. The filtered schemes 
uses the accurate method at more locations than the hybrid discretization. This 
occurs despite the fact that the hybrid scheme used known regularity results to 
choose the locations of the more accurate method, whereas the filtered scheme 
determines the discretization from the equation. 

Computation speed. Table 2 presents the number of Newton iterations, compu- 
tation time, and maximum error in each of the computed solutions for the hybrid 
and filtered schemes on a 17 point stencil. As a benchmark, the error in the stan- 
dard scheme is also presented; it has been shown previously that this scheme is 
much slower than the Newton solver for the hybrid discretization [FOllb]. Even 
on this relatively narrow stencil, the accuracy of the filtered scheme is close to — 
and in some cases better than — the accuracy achieved using the slow, formally 
second-order standard scheme. Overall, there is no appreciable difference in accu- 
racy between the results obtained using the hybrid and the filtered scheme, though 
the filtered scheme is slightly more accurate in most of the examples. (With our 
particular choice of the parameter e{h, dO) it is slightly less accurate on the blow-up 
example. Accuracy can be improved by allowing this parameter to increase). Of 
course, the filtered scheme has one big advantage over the hybrid scheme because 
it comes with a convergence proof, which gives us confidence that the method will 
continue to perform correctly in other examples. 
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(B) 

Figure 6. Illustration of the weights in the filtered scheme (6a), 
and the hybrid scheme (6b). The ordering is the same as in Fig- 
ure 4. In the second example top figure, green indicates a weight 
of 5% on the monotone scheme. In the third and fourth figures, 
green indicates the monotone scheme, while yellow indicates the 
accurate scheme. 

5. Conclusions 

We constructed and implemented a convergent, higher order accurate scheme for 
the Monge- Ampere equation. We extended the convergence theory of Barles and 
Souganidis by considering the more general class of nearly monotone schemes. This 
new convergence proof applies in general to the class of nonlinear elliptic PDEs, 
and requires only an elliptic scheme as a foundation. 

The combined schemes are called filtered finite difference approximations. The 
filtered scheme chooses between a convergent elliptic scheme and a more accurate 
scheme. The selection principle is based on filtering the difference between the 
elliptic scheme and the more accurate scheme, reducing to the elliptic scheme when 
the difference is large, and the accurate scheme when the difference is small. 

The theory ensures, and computational results verify, that solutions of this 
scheme converge to the viscosity solution of the Monge- Ampere equation even in 
the singular setting. 

The accuracy of the filtered schemes was as good as the accuracy of the other 
methods. Formal 0{h'^) accuracy was attained on the smooth example, while on 
singular examples the accuracy decreased as expected. Newton's method resulted in 
a fast solver, with the same number of iterations (2-11, depending on the solution) 
as for the monotone method and the hybrid method. This is much faster than other 
types of solvers, which typically break down on singular solutions, or require more 
iterations on larger problems. 
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Max Error, Example 
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9 Point 
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Monotone 
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1.04 X 10- 
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■5 
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0.26 X 10- 


-5 
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-5 
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0.26 X 10- 


-5 
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■5 
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-5 
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0.06 X 10" 


■5 
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3.48 X 10-5 
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■6 


1.08 X 10-5 
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Max Error, Example 
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Filtered 
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4 
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-4 
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-4 
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0.46 X 10" 
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-4 
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21.53 X 10-* 


0.23 X 10" 


-4 
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0.31 X 10" 


-4 


2.37 X 10-* 


0.28 X 10" 


-4 



Max Error, Example with Blow-up 
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1.74 X 10" 
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1.74 X 10" 


-3 


1.74 X 10- 


3 


1.74 X 10" 


-3 


1.74 X 10" 


-3 


63 


0.86 X 10-3 


0.59 X 10- 


3 


0.59 X 10- 


-3 


0.59 X 10- 


■3 


0.59 X 10- 


-3 


0.59 X 10- 


-3 


127 


0.83 X 10-3 


0.20 X 10" 


3 


0.35 X 10" 


-3 


0.20 X 10" 


■3 


0.20 X 10" 


-3 


0.20 X 10" 


-3 


255 


0.83 X 10-3 


0.15 X 10- 


-3 


0.33 X 10- 


-3 


0.13 X 10- 


3 


0.16 X 10- 


-3 


0.08 X 10- 


3 


361 


0.83 X 10-3 


0.17 X 10- 


3 


0.33 X 10- 


-3 


0.13 X 10- 


3 


0.15 X 10- 


-3 


0.08 X 10- 


3 



Max Error, C^'^ (Lipschitz) Example 
N 9 Point 17 Point 33 Point 





Monotone 


Filtered 
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Filtered 




Monotone 


Filtered 




31 


11.83 X 10-3 


10.42 X lo- 


-3 


3.56 X 10-3 


4.16 X 10- 


-3 


1.61 X 10-3 


1.23 X 10- 


-3 


63 


11.10 X 10-3 


ll. 56 X 10- 


-3 


3.49 X 10-3 


2.82 X 10- 


-3 


1.65 X 10-3 


1.90 X 10- 


-3 


127 


11.80 X 10-3 


10.97 X 10- 


-3 


3.45 X 10-3 


2.83 X 10- 


-3 


1.64 X 10-3 


1.11 X 10- 


3 


255 


10.47 X 10-3 


11.03 X 10- 


-3 


3.46 X 10-3 


3.06 X 10- 


-3 


1.64 X 10-3 


1.15 X 10- 


■3 


361 


10.40 X 10-3 


10.37 X 10- 


-3 


3.45 X 10-3 


3.22 X 10- 


-3 


1.64 X 10-3 


1.12 X 10- 


■3 



Table 1 . Accuracy of the monotone and filtered schemes. 



For the convergence theory, we require d9 — >■ 0, but in practice, using the filtered 
scheme means the accuracy is not affected by d0, so we can use the narrow stencil 
scheme and still obtain accuracy corresponding to the regularity of the solution 
(e.g. 0{h?) on smooth solutions, 0{h) on moderately singular solutions). The only 
exception was the most singular example, which was not a viscosity solution, in any 
case. 

In summary, the new method combines all the advantages of our previous meth- 
ods: speed of solution, accuracy, and a proof of convergence, while also allowing 
the use of a narrow stencil scheme in practice. 
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12 


101.5 


108.7 


0.04 X 10- 


-3 


0.13 X 10"^ 


5.42 X 10-3 










Lipschitz Example 










Iterations 


CPU Time (s) 






Max Error 




N 


Hybrid 


Filtered 


Hybrid 


Filtered 


Hybrid 




Filtered 


Standard 


31 


6 


7 


0.2 


0.2 


3.57 X 10" 


-3 


4.16 X 10-3 


5.19 X 10-3 


63 


6 


9 


1.0 


1.4 


3.49 X 10- 


-3 


2.82 X 10-3 


2.86 X 10-3 


127 


8 


8 


6.6 


6.4 


3.45 X 10- 


-3 


2.83 X 10-3 


1.54 X 10-3 


255 


9 


10 


47.1 


38.2 


3.46 X 10- 


-3 


3.06 X 10-3 




361 


10 


9 


155.8 


81.7 


3.45 X 10" 


-3 


3.22 X 10-3 





Table 2. Newton iterations, computation time, and accuracy for 
the hybrid and filtered schemes on a 17 point stencil. 
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